A simple dynamical model with history dependence for a sand-pile experiment 
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A lattice dynamics model is proposed for the history de- 
pendence observed in sandpile experiments. The dependence 
of the stress distribution on the preparation of the sandpile is 
explained as a dependence of certain attractors on the prepa- 
ration of the system. The model has three phases, but the 
history dependence is shown to exist only in the phase where 
a perturbation is amplified selectively rather than globally 
when propagating in the downflow direction. The condition 
for this history dependence is given in terms of the spatial 
Lyapunov exponent. 

PACS number(s): 45.70.Cc, 05.45.-a, 05.45.Ra 



Physical systems exhibiting a "history dependence" in 
the sense that their state is not solely determined by ex- 
isting environmental conditions but dependent on how 
the system was externally driven in past are not uncom- 
mon. Examples of such systems can be found among 
glasses, polymers, and granular materials, and of course, 
in biological systems. Although hysteresis between two 
states is the simplest form of history dependence, here we 
examine cases where the evolution of a system may lead 
to many different states. Consequently, some variable(s) 
for the internal state can be considered a type of 'mem- 
ory'. In a dynamical system, each memory state can be 
interpreted as a different attractor. Hence it is necessary 
to study how each attractor is selected and how feasi- 
ble it is for one attractor to switch to another attractor. 
In order to do so, we take a spatially extended system, 
and discuss the macroscopic features of the history de- 
pendence. 

A simple example of history dependence is given by 
recent sand-pile experiments In these experiments, 
it is found that the pressure profile formed by the stress 
chains in a sand pile strongly depend on how the sand 
beads are piled. In the present Letter, we introduce a 
simple toy model inspired by the experiment, and dis- 
cuss some characteristic features of a system with history 
dependence. 

In our toy model we take a mesoscopic approach by us- 
ing a coarse-grained stress field instead of taking a vari- 
able for each bead position. A 'stress variable' is 
assigned to each site located on a two-dimensional trian- 
gular lattice Note that due to gravity, the stress 



in a sand pile increases from upper to lower layers. Each 
site is associated with an intrinsic gravitational weight 
such that, when excluding the stress from higher layers, 
the stress variable relaxes according to the weight of the 
site given by /3y. The stress of a layer is transferred 
downwards by distributing the stress on each of its sites 
to the right and left neighbor sites in the next layer. The 
amount of stress transferred from the upper sites that can 
be supported by each of the two lower sites is assumed 
to depend on the existing stress values of the lower sites. 
The larger the existing stress value of a site is, the more 
effect the stress from the upper site will have. 
Based on these arguments, we introduce the following 
model of coupled differential equations: 



_ dxjj(i) 
~dt 



-Xij(t) + Pij + X hj -i(t)- 



i(t)Y 



+ X h+ ij-i(t)- 



with L = 



Xij {t)] a + [x mj (t)} a 

i for even j 
i — 1 for odd j 



[a*-ii(*)] a 
(1) 



where a is a positive parameter that characterizes the 
tendency of a site with a higher existing stress value to 
be affected more strongly by an upper layer. Here the site 
(Ji,j — 1) is the upper left of the site and the site 
is the upper right. The time scale r can be set 
at 1, by rescaling the time. We consider the homogeneous 
case here with a constant /Jy = (3 (which can be set to 
1 by rescaling the variable x). The lattice size, unless 
mentioned otherwise, is chosen to be N for the horizontal 
direction with periodic boundary conditions, and M for 
the vertical direction (sites in both the horizontal and 
vertical directions are numbered from and hence j = 
M — 1 is the bottom layer || ) . 

In our model the 'weight' of the site (i,j — 1) is sup- 
ported by the two neighboring sites in the jth layer, 



with the ratio of xf j /(xf j + 



-ki 



). Corresponding to 



the 'conservation law' of momentum flux, the relation 

= {xj-i)i - {xj)i + {J3j-i)i holds, where ( ) 4 de- 
notes the average over the horizontal direction for the jth 
layer, (xj)i = j?J2i x v and (&-!>* = wE*/%-i = 1- 



\Xj)i — N x i 
Since the equation is governed by the relaxation term, 
the stress variables Xij(t) are attracted to a fixed point 
solution x*j . As will be shown, there is a huge number of 
stable fixed point solutions for the steady state equation 
if a > 0.7. For a fixed point solution, the conservation 



law implies that {x*)i = (af?_i)i + (/3j_i)j 
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Now we study the history dependence in our model by 
recalling the sandpile experiments. For a sand pile, it is 
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known |lj that when the sand is piled from a small hole, 
the pressure distribution at the bottom has a minimum in 
the center, while it has a maximum at the center when the 
sand is scattered during the piling. With this experiment 
in mind, we choose the following operations to prepare a 
'sand pile'. 

In our model, we represent the process of piling as the 
addition of a weight fyj = 1 at a site, while the stress 
variable Xij then evolves from 0. Before piling a bead 
at a site, 0ij = and iy = were assigned. Then, in 
order to simulate a piling process from bottom to top, at 
each step, we choose a horizontal site i randomly from a 
Gaussian distribution with the standard deviation W x 
(N/2) and took the largest value j such that the site 
is not vacant. If both the neighboring lower sites 
are occupied, the weight f3 is added at that site. If at 
least one of the sites of the neighboring lower sites is 
vacant, this (meso-scopic) bead falls down to the lower 
vacant site. By adding weights successively, the dynamics 
is iterated to lead to a fixed stress field. 

The case with W <C 1 corresponds to a piling process 
from a thin hole, while on the other hand the case with 
W > 1 corresponds to a piling process where the beads 
are scattered. 

Figure [l] shows the results of the stress field for (a)W = 
0.1, and (b)VK=10. Note first, stress accumulated 
on a few lattice sites forming a downflow stream. These 
sites correspond to the stress chain observed in the sand 
pile experiments. The accumulation of the stress itself is 
obvious, since our model has a tendency of more stress 
at a site with a larger stress value. In fact, the formation 
of stress chains is commonly observed for a > 0.7. 

As shown in the figure the stress chain is distributed 
around the edge for the case (a) (Figure 0(a)) , while in 
the case (b), it is distributed around the center. The 
stress distribution at the bottom plate is plotted in 
Fig.[l](c), obtained from the average of 100 samples. It has 
a dip at the center for the case (a), while it is unimodal 
for (b). The stress distribution with the dip is observed 
independent of the system size and fits a scale-invariant 
form .These results agree well with the observations in 
sand-pile experiments Q and particle simulations. Note 
that the history-dependence shown here is generally ob- 
served for a > 1. 

As a second demonstration of history dependence, we 
study how the stress pattern depends on the time scale 
of preparing a sand pile. Differences in time scales are 
introduced as follows: first the system is again prepared 
in a "no bead state ", i,e„ a state with = x^ = 
for all sites. Then, every T steps ( untill the top layer 
j = is reached), we set f3 = 1 for all the sites in the 
lowest vacant layer simultaneously. The stress field x^ 
of a new site is set to a random value in the range [0, 0.1] 
to represent small fluctuations. Thus, with the help of 
eq.(Q), a fixed stress pattern x*j is obtained. 

Through this process, stress chains are again formed. 
We have computed the number of stress chains ((N )) by 
counting the sites that satisfy x*j > {x*)i = j + 1 in the 



bottom layer, where ((..)) is the ensemble average over 
500 trials. When the piling process is fast (T < 10), the 
number of stress chains is small, and in fact the stress 
is accumulated at few sites in a lower layer. For larger 
T the number of stress chains is larger, and the weight 
is supported by more sites |(|. This dependence on the 
time scale is commonly observed for a>l. The time 
scale T w 10 is related to the intrinsic time scale for 
the relaxation of the coupled system, as will be discussed 
later. 

Now we discuss the origin of this history dependence in 
terms of dynamical systems theory. First, the coexistence 
of multiple attractors is confirmed by choosing a variety 
of initial conditions with x^j <E [0,1]. Then after the 
relaxation is completed we compute the patterns of the 
fixed point attractors. Depending on the value of a, we 
have the following three phases. 

(T)a < 0.7: The system is always attracted to a sta- 
ble homogeneous state. A single stable fixed-point with 
x*j = j + 1 (independent of i) exists. 

(II) 0.7 < a < 1: Besides the homogeneous solution, 
many stable fixed point patterns appear (Fig.^|(a)). In 
the upper region, the patterns are almost homogeneous 
while in the lower region they are striped (zigzag-like) 
with some dislocations. For random initial conditions, 
the homogeneous state is not realized unless we take an 
initial condition in its vicinity. 

(III) a > 1: The homogeneous solution becomes un- 
stable at a = 1 + + 1), and further downflow it is 
always unstable. When following the layers down, the 
behavior changes from an almost homogenous state, to a 
stripe pattern as in phase II, and then to patterns where 
stresses are accumulated into fewer sites forming inho- 
mogeneous patterns as shown in Fig.||(b). 

Although both phases II and III have multiple fixed- 
points, the history dependence discussed so far is ob- 
served only for phase III. One important difference be- 
tween phases II and III lies in the response of each at- 
tractor to a perturbation applied at the top layer. 

In phase II, there is a switch from one attractor to 
another when applying a tiny perturbation for at least 
0(1) time units. A local perturbation in an upper layer 
spreads out horizontally in lower layers. Fig.|^(a) shows 
x^ and the variation Ax|- caused by the perturbation. 
As can be seen, x^ are altered for almost all sites. 

In phase III, a perturbation needs to be applied for 
some time in order to cause a switch to a different at- 
tractor. Furthermore, the length with which the pertur- 
bation needs to be applied increases as the perturbation 
amplitude becomes smaller. Here, a perturbation is not 
transferred smoothly to all sites, but is localized around 
a few stress chains. As shown in Fig.[| the sites whose 
values are altered are located around the stress chains, 
especially in the case of lower layers. The change of x*j 
is temporally intermittent. After a perturbation is ap- 
plied, the stress values stay almost constant but then for 
some sites they show a sudden increase, corresponding to 



a rearrangement of the stress chains. 

In order to study the relationship between the response 
to a perturbation and the history dependence, we carried 
out the following simulations: After the system reaches 
a fixed point, we apply a perturbation at the top (0-th) 
layer to change /3q = (1, 1, • • ■) to 0o + 5f3 over r steps, 
where 5(3q is chosen randomly over [— e,e] and fixed. Af- 
ter waiting for the system to settle down to a fixed point 
attractor, we checked whether it is identical to the orig- 
inal attractor or not. By taking 100 samples, the ratio 
of attractor switching is computed as a function of the 
perturbation strength e and the perturbation duration t. 
This is shown in Fig.|| for phase III. As the duration time 
t for the application of the perturbation increases, its ef- 
fect accumulates causing a switch from one attractor to 
another. The required time duration for a given switch- 
ratio roughly scales as 1/e. Hence the history dependence 
in Fig.^j is a natural consequence of the above time scale 
dependence. Note that for phase II, a dependence on r 
or e is not observed since almost all perturbations kick 
the system to a different attractor. 

The response of a system against tiny perturbations is 
generally studied by means of a Lyapunov analysis. Or- 
dinary Lyapunov analysis, however, is not useful here, 
since all the fixed-points studied so far are linearly sta- 
ble. Rather, in an open-flow system with a unidirectional 
interaction in space, the co-moving Lyapunov exponent 
measuring how a perturbation is amplified along the spa- 
tial direction is important Q . 

In order to discuss the stability along the spatial di- 
rection, we first obtain the relationship between the 
fixed points at the (j — l)-th layer and the j-th layer 
x*j = f(x*j-x,x*j) by equating the r.h.s. of eq.(Q) with 
0. This gives a spatial map from x* j_i to x*j. The am- 
plification rate of the perturbation 5x*_ 1 at the (j — 1)- 
th layer is given by 8x* = Mj5x*_ l , with Mj the Jacobi 
matrix for the spatial map. 

The rate at which a perturbation expands in the down- 
flow direction is then given by the maximal Lyapunov 
exponent of the spatial map as X^ p — log(/ij)/j, where 

/ij is square root of the maximal eigenvalue of Q^Qj, 
with Qj — MjMj-i ■ ■ ■ Mi'i and j starting from the 0-th 
layer. The local expansion rate from one layer to the 
next is computed as the logarithm of the square root of 
the maximal eigenvalue of MjMj. We have plotted A* p 
in Fig.|| as a function of the layer, and X l ° c in its inset. 

In phases II and III (a = 0.8, 1.2) A^ oc - is positive im- 
plying that a perturbation in one layer is amplified to 
the next layer. The behavior of A* p , however, is differ- 
ent between the two phases. In phase II (for a = 0.8) 
A* p increases with the layer j and in fact the attractors 
of phase II are convectively unstable Q. On the other 
hand, in phase III (for a— 1.2), X s ? first increases sharply 
and then decreases towards as j is increased. Hence, a 
perturbation in the top layer may not reach the bottom 
layer if the number of layers is sufficiently large. 



This property is important for explaining the history 
dependence shown in the piling process. Since the local 
exponent is positive, a perturbation in the piling pro- 
cess is amplified to the next layer, and thus leads to a 
variety of stress patterns. As the sand is piled, the dis- 
tance between the top and the bottom layers increases, 
and a perturbation by the addition of a layer at the top 
has less influence on the bottom layers because A* p ap- 
proaches zero. The lower layers are thus 'protected' from 
the changes at the top layers. Hence a given stress pat- 
tern selected in the initial stage is 'memorized'. This 
leads to the history dependence in phase III, while for 
phase II, a perturbation in the top layer is continuously 
amplified to the bottom (since lim^oo A^ p > 0) and there 
is no history dependence. 

In conclusion, we have constructed a simple model that 
describes the history dependence of the pressure distri- 
bution in a sand-pile experiment. A condition for the 
history dependence is discussed from the viewpoint of a 
response against perturbations. When there is a strong 
constraint on the attractors that can be reached by a 
perturbation to a given attractor, history dependence 
is observed, while if the perturbation influences all sites 
globally, history dependence is not possible || . This con- 
dition for history dependence was further characterized 
by analyzing the spatial Lyapunov exponent. 

The history dependence we have observed should be a 
general property of systems with the self-reinforcement 
of stress (as represented by the x a term in eq.(|l|), but 
several other forms may also give this type of history 
dependence). How such a reinforcement term appears 
from a microscopic dynamics of granular matter is left for 
future work, while it will also be important to generalize 
the present study to other history dependent phenomena 
including biological systems. 
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(a)W=0.1 V (b)W=10.0 




FIG. 1. Sand-pile simulations with 

(a)W = 0.1(b)W = 10.0. White pixels indicate large xj^i.e. 
x*j > (x*)i). (c) Distribution of x *m-i obtained by averaging 
100 samples for W = 0.1 (red line) and W=10.0(green line). 
The data are for N = 128, but the results for N = 256 are 
qualitatively identical after scaling. 
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FIG. 2. The ratio of the number of stress chains and N 
plotted against the time scale T of the piling process described 
in the text. The vertical lattice size is 120, and the horizon- 
tal size is 64(d), 128(o),200(«). The results were obtained by 
averaging 500 runs. 
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FIG. 3. The left columns in (a) and (b) show the stress vari- 
ables x*j plotted in gray scale for a system size of 32 x 200. 
(a) a = 0.9 and (b)a = 1.2. The right columns depict 
Ax*j — \x*j —x*j | in gray scale, where the x*j is the fixed point 
reached after a perturbation is added as /3o = (1,1,- ■ •)— >/3o+5/3, 
with 5/3 G [-0.00001,0.00001]. 
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FIG. 4. The fraction of attractor switches due to a per- 
turbation with an amplitude e £ [0.0001,0.1]. and a time 
duration r in the top layer. The lattice size is 32 x 200 and 
a = 1.2. 
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FIG. 5. Spatial Lyapunov exponent A° p and local spatial 
exponent (inset) for a = 0.8(x) and 1.2(o). The lattice size 
is 32 x 200. 



